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ABSTRACT 

The physical properties inferred from the spectral energy distributions of z > 3 galaxies have been 
influential in shaping our understanding of early galaxy formation and the role galaxies may play 
in cosmic reionization. Of particular importance is the stellar mass density at early times which 
represents the integral of earlier star formation. An important puzzle arising from the measurements 
so far reported is that the specific star formation rates (sSFR) evolve far less rapidly than expected in 
most theoretical models. Yet the observations underpinning these results remain very uncertain, owing 
in part to the possible contamination of rest-optical broadband light from strong nebular emission lines. 
To quantify the contribution of nebular emission to broad-band fluxes, we investigate the spectral 
energy distributions of 92 spectroscopically-confirmed galaxies in the redshift range 3.8 < z < 5.0 
chosen because the Ha line lies within the Spitzer/IRAC 3.6/im filter. We demonstrate that the 
3.6/im flux is systematically in excess of that expected from stellar continuum alone, which we derive 
by fitting the SED with population synthesis models. No such excess is seen in a control sample of 
spectroscopically-confirmed galaxies with 3.1 < z < 3.6 in which there is no nebular contamination 
in the IRAC filters. From the distribution of our 3.6/mi flux excesses, we derive an Ha equivalent 
width distribution and consider the implications both for the derived stellar masses and the sSFR 
evolution. The mean rest-frame Ha equivalent width we infer at 3.8 < z < 5.0 (270 A) indicates that 
nebular emission contributes at least 30% of the 3.6/im flux and, by implication, nebular emission 
is likely to have a much greater impact for galaxies with z ~ 6 — 7 where both warm IRAC filters 
are contaminated. Via our empirically-derived equivalent width distribution we correct the available 
stellar mass densities and show that the sSFR evolves more rapidly at z > 4 than previously thought, 
supporting up to a 5x increase between z ~ 2 and 7. Such a trend is much closer to theoretical 
expectations. Given our findings, we discuss the prospects for verifying quantitatively the nebular 
emission line strengths prior to the launch of the James Webb Space Telescope. 

Subject headings: galaxies: formation - galaxies: evolution - galaxies: starburst - galaxies: high 
redshift - ultraviolet: galaxies - surveys 



1. INTRODUCTION 

Through detailed photometry of Lyman break galaxies 
(LBGs) undertaken with the Hubble and Spitzer Space 
Telescopes, much has been learned regarding the physi- 
cal properties of galaxies beyond a redshift z ~3. Stellar 
masses and star formation rates have now been inferred 
from broadband photometric spectral energy distribu- 
tions (SEDs) for thousand s of galaxies spanning the red- 
shift r ange 3 < z < 7 (e.g ., lEgami et al. 2005; Evlc s et al 
2005 1 : lLabbe et all l2006t lEvles et all l2007t IStark et al 
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Labbe et al.l l2010albl 

20121 : iReddv et al.ll2012bb iCurtis-Lake et al.ll2012f l The 
stellar mass density derived from these studies has proven 
a useful integrated constraint on the contribution o f 
galaxies to reionization (e.g., iRobertson et al.l l2010h . 
while the evolution of physical properties has provided 
insight into the processes which govern the assembly of 
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early galaxies (e.g.. lFinlator et aTll2011l : lDave et al.ll201lL 
[2012TI . 

A potentially significant puzzle has recently emerged 
from these studies through measurement of the specific 
star formation rate (sSFR) at z > 2. Current observa- 
tions demonstrate that between z ~ 2 and z ~ 7, the 
sSFR in galax ies of fixed stellar mass does not evolv e 
strongly (e.g., IStark et al.l l2009t IGonzalez et al.llMof ). 
with recent estimates indicating at m ost a factor of 
two increase between z ~ 2 and 7 (e.g., iBouwens et al.l 
l2012bt IReddv et al1l2012bh . This is in contrast to sim- 
ple expectations fro m semi-analytic models and numer- 
ical s i mulat ions (e.g.. IWeinmann et a,l.ll201 it iDave et al.l 
120111 12012ft which predict that the sSFR should closely 
match the inflow rate of baryonic material. As this 
mass inflow rate is thought to increase with redshift as 
M/M ~ (1 + z) 2 - 25 (jNeistein fc Dekelll200l iDekel et al.1 
2009), we should expect nearly a 10 x increase in sSFR in 
galaxies of fixed stellar mass over 2 < z < 7, in marked 
contrast to the observations. 

The physical cause of the discrepancy associated 
with the sSFR evolu tion remains unclear. As dis- 
cussed previously (e.g.,[ Bouchc et al. 2010; Dutton ~et al.1 
20101: IWeinmann et all 120111 : IDave et all 1201 lL 120121 : 



Reddv et all l2012b|) , a plateau in the redshift depen- 



dence of the sSFR would suggest star formation is more 



2 



Stark et al. 



inefficient at z > 6 than at z ~ 2. Various physical 
processes might be invoked to impede star formation, 
such as the inefficient formation of molecular hydrogen 
in low metallicity galaxies (e.g., iRobertson fe Kravtsovl 
2008: lGnedin et al. 2009: K rumholz fe De kcl 2012) or an 
increase in the efficiency with which cold gas is removed 
via large-scale outflows. 

Irrespective of any mechanism that might inhibit star 
formation at early times, it is difficult to reconcile such 
an inefficiency with the notion that early galaxies pro- 
yide the ionizing pho tons responsible for reionization 
((Robertson et alll2010|) . For example, the steep faint end 
slope of the ultraviolet lum inosity function (UV LF) at 
z > 6 (jBouwens et al.ll2011[ ) implies that star formation 
in low mass dark matter halos becom es more efficient at 
earlier times (e.g.. iTrenti et al1 l2010). in contrast to the 
implications of the sSFR measurements. 

Given these difficulties, it is prudent that we recon- 
sider the accuracy of the data that is used to infer the 
sSFR and its evolution. The two basic ingredients are 
the star formation rates and the stellar masses. The 
z > 4 measurements have indeed c hanged since the orig i- 
nal articles (e.g. JStark et al.ll2009l: [Gonzalez et al.ll2010h . 
mostly as a result of improved dust corre ctions following 
improv ed near-infrared photometry (e.g.. lBouwens et al.l 
20123). The new dust corrections have served to in- 
crease the z ~ 4 sSFR measurements by a factor ~ 
2. However, since negligible extinction is inferred at 
z ~ 6 — 7, the sSFR sti ll remains constant over 4 < z < 7 
(|Bouwens et alll2012bf) although a factor of 2 higher than 
at z ~ 2 - 3 (|Reddv et al.l 1201 2bD . 

A potentially more important problem is the possible 
contribution of rest-frame optical nebular emission lines 
(e.g., [O II], [O III], Ha) to the broad-band fluxes used 
to infer the stellar masses. Such emission lines, could 
significantly affect the inferred amplitude of a Balmcr 
Break, leading to an overestimate of the stellar mass and 
thereby an underestimate of the sSFR. Figure 1 illus- 
trates how the various nebular emission lines contami- 
nate the key photometric filters as a function of redshift. 
Beyond z ~ 4, the key filters of interest in the determina- 
tion of stellar masses are the Spitzer/IRAC warm bands 
at 3.6/zm and 4.5/im. It is particularly striking that, 
at z > 5, the strongest rest-frame optical nebular lines 
([OIII]A5007 and Ha) contaminate both Spitzer/IRAC 
filters. Although many z > 5 galaxies are detected 
with Spitze r (e.g. lEgami et al l 120051: lEvles et "ail 120051: 
Stark et al.l 120091: lLabbe et all boiM bl: iGonzalez et al.1 
20101 I2011al : iRichard et ail 1201 lbf l. contamination by 
nebular emission could significantly affect the interpre- 
tation of their SEDs. 

Accounting for nebular emission in the SEDs of high 
redshift gal axies has been considered by several earlier 
work s (e.g.JSchaerer fc: de Barrosl l2009l (20101: lOno et all 
120101: Ide Barros et al.l 120121 " In general terms, their ap- 
proach has been to use 'forward modeling' techniques 
based on adding nebular emission contributions to stel- 
lar population synthesis models in order to demonstrate 
the possible implications of its inclusion. However, such 
"nebular+stellar" model fits cannot provide a precise un- 
ambiguous measure of nebular contamination for several 
reasons. Firstly, there are numerous uncertainties in how 
the contribution of nebular emission should be added. 



These include the nebular extinction law and ionizing 
photon escape fraction. Secondly, for galaxies at z > 5, 
for which there is no uncontaminated measure of the stel- 
lar continuum (Figure 1), the uncertainties are particu- 
larly large. Finally, and perhaps most importantly, with- 
out a spectroscopic redshift, addressing both the nebular 
contamination and the photometric redshift of the galaxy 
from the same photometric data leads to great uncertain- 
ties; there is no a priori indication of which photometric 
bands are contaminated by nebular emission. 

Fortunately, by virtue of our deep Keck spectroscopic 
survey (Stark et al. 2010, Stark et al. 2011, Jones et 
al. 2012, Schenker et al. 2012) and our nebular+stellar 
population synthesis code (Robertson et al. 2010), we 
can use the availability of HST-Spitzer SEDs to make 
progress in addressing this issue. While the question 
of contamination by nebular emission at z > 5 must 
await the infrared spectroscopic capabilities of James 
Webb Space Telescope, we can test our spectroscopic 
range 3.8 < z < 5.0 for contamination by Ha in the 
Spitzer/IRA C 3.6/nrn broadban d filter. Our approach fol- 
lows that of IShim et all (|2011[) who demonstrated that 
galaxies in this redshift window are typically signifi- 
cantly brighter at 3.6/im than at 4.5/zm. By comparing 
their flux density at 3.6/mi to that e xpected from stellar 
continuum alone, IShim et"ai1 (|2011l ) argued that many 
galaxies at z > 4 show evidence for strong Ha emission, 
with typical EWs significantly greater than those seen at 
z~2. 

Here we seek to apply a similar techniq ue to our spec- 
troscopic sample ([Stark et al.ll2010l 1201 ll ) with the goal 
of estimating the distribution of Ha equivalent widths 
present in galaxies at 3.8 < z < 5.0. Equipped with this 
external constraint on the strength of nebular emission, 
we can then determine how stellar masses and the sSFR 
of z > 4 galaxies are likely to be altered by emission line 
contamination. In particular, we will explore whether 
our estimated degree of nebular contamination could be 
sufficient at the highest redshifts to permit a rapid rise 
in the redshift-dependent sSFR as expected from theo- 
retical models. 

The present paper is organized as follows. In §2, we 
discuss the selection of the spectroscopic sample used in 
our analysis. In §3, we introduce the details of our SED 
fitting procedure used to estimate the strength of nebu- 
lar emission lines in the various filters. In §4 we use our 
spectroscopic sample to estimate the equivalent width 
distribution of Ha in the redshift range 3.8 < z < 5.0 
and then use these measurements to assess the impact of 
nebular emission on the derived stellar masses and star 
formation rates of z > 4 galaxies. In §5, we discuss the 
impact that our findings have for the evolution in the in- 
tegrated stellar mass density and specific star formation 
rates of galaxies at z > 4. 

Throughout this paper we adopt a A-dominated, flat 
universe with = 0.7, flju = 0.3 and Hq = 

70 /^okms" 1 Mpc -1 . All magn itudes are quoted in the 
AB system (jOke fc Gunnlll983H 

2. DATA 

In this paper, we will focus our analysis on the in- 
terpretation of broad-band spectral energy distributions 
for a z > 3 sample with known spectroscopic redshifts. 
The spectroscopic sample is drawn from earlier papers 
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Fig. 1. — Emission line contamination of broadband photometry. 
Colored stripes denote redshift ranges over which emission lines 
contaminate the K s -band (dark blue), IRAC 3.6 fim (yellow), and 
IRAC 4.5 fim (red). Ho emission is expected in the 3.6/^m filter 
at 3.8 < z < 5.0. Note that at 5 < z < 7, both IRAC filters used 
for measuring stellar masses are contaminated by emission lines. 
Beyond z ~ 7, only the 4.5/xm filter is contaminated by strong 
nebular emission. 



(jStark et al.l 120101 l201ll) . Full details can be found in 
these articles, but we offer the reader a brief summary 
here. Spectroscopy of z > 3 LBGs in the two GOODS 
fields was undertaken at t he Keck Observat ory using 
the DEI MOS spectro graph dFaber et al1 l2003t). As dis- 
cussed in iStark et al.l (|2010l ). LB Gs were selected using 
standard 'dropou t' criteria (e.g., iBouwens et al.l l2007t 
IStark et al.l [2009) to a limiting magnitude of 2:350 — 27 
using the GOODS v2 public photometric catalogs (e.g., 
iGiavalisco et al.1 l2004f ) . Taking advantage of a similar 
spectroscopic campa ign undertaken using the FORS2 in 
GOODS-South (e.g- lVanzella et al.l kOOQ). we retrospec- 
tively constructed a VLT sample using the same pho- 
tometric criteria. The combined Keck plus VLT survey 
comprises 157 galaxies in the redshift range 3.8 < z < 5.0 
which satisfy the dropout criteria. As we will discuss be- 
low, only a subset of these will be used in our analysis. 

A key requirement for the derivation of stellar masses 
and specific star formation rate is precise broad band 
photometry from which SEDs for galaxies of known 
spectroscopic redshifts can be determined. In GOODS- 
South, we use the public release of the Wide Field Cam- 
era 3 (WFC3) ima ging from the CANDELS Multi-Cycle 
Treas ury Program (|Grogin et al.ll201ll:lKoekemoer et al.l 
120111) and our own reduction (see lMcLure et alll2011f) of 
the Early Release Science (ERS) campaign (e.g., Wind- 
horst et al. 2011). Colors were computed with respect 
to the zg5o flux using matched apertures with up-to-date 
zero points, and total WFC3 magnitudes were derived by 
combining the measured colours with the total zsso-band 
flux. Ks- band photometry is taken from deep ISAAC 
imaging (iRetzlaff et al.l [2010h following the procedure 
discussed in IStark et al.l (|2009l ). For GOODS-N, we use 
near-infrared imag ing obtained from CFHT/WIRCAM 
(|Wang et al.ll201oh . 

The rest-frame optical at z > 4 is probed by the deep 
Spitzer/IRAC (Fazio et al. 2004) imaging of GOODS- 
S and GOODS-N (Dickinson et al. 2012, in prep). In 
particular, the 3.6/mi (hereafter [3.6]) and 4.5/im (here- 
after [4.5]) are the most useful, as the longer wavelength 
filters are typically not sensitive enough to detect most 
z > 4 galaxies. As in Stark et al (2009), we focus primar- 
ily on the subset of ACS-selected galaxies whose IRAC 



fluxes are not contaminated significantly by neighboring 
sources. The IRAC magnitudes are measured in aper- 
tures 2.4 arcsec in diameter and to account for flux falling 
outside this aperture, we apply a 0.7 mag aperture cor- 
rection derived from a sample of isolated point sources. 
Recognizing that selecting only isolated IRAC sources 
limits the size of our eventual sample, we included IRAC 
flux measurements fo r galaxies in GOODS-South from 
the MUSIC catalog (jGrazian et al.l 120061 : ISantini et all 
2009). These fluxes rely on a deconfusion procedure to 
extract fluxes from sources with contaminating neigh- 
bors. A comparison between the two photometry meth- 
ods reveals consistency for our isolated sample, with a 
standard deviation of 0.19 mags and no systematic off- 
set. 

In total, we have 92 galaxies in the range 3.8 < z < 5.0 
with measured IRAC photometry. In our analysis, we 
will focus on the subset of 45 galaxies with confident 
(>5<t) 4.5/im detections (see Table 1), as without an ac- 
curate measure of the 4.5/im flux it is impossible to infer 
the expected stellar continuum from population synthesis 
models. The objects with deconfused GOODS MUSIC 
photometry make up 60% of the final sample. 

Some caution must be exercised when applying infer- 
ences from a spectroscopic sample to the parent photo- 
metric population. In attempting to infer the typical 
level of rest-optical nebular contamination, we must be 
particularly careful that we do not bias our sample to- 
ward strong Lya emitting galaxies, a population which 
might have larger than average sSFR and Ha EW. For 
faint galaxies (zgso > 25.5), the spectroscopic sample of 
Stark et al. (2010) is indeed biased toward Lya emit- 
ters. But by requiring a 5cr detection in the [4.5] band, 
we limit our sample to brighter systems (average zsso- 
band magnitude of 25.0) for which we are more complete 
spectroscopically. Indeed the percentage of galaxies for 
which we measure Lya in emission is actually only 46%, 
highlighting the fact that many galaxies in this bright 
subset are instead confirmed via the combination of Lya 
in absorption and metal absorption lines (e.g., Jones et 
al. 2012). Furthermore, the fraction of galaxies in this 
subset with strong (EW > 50 A) Lya emission (5%) is 
similar to that measured for galaxies in this Myy and 
redshift range (6% in Stark et al. 2010) adding confi- 
dence that the sample we use in this paper is not likely to 
be strongly biased toward nebular emitters and appears 
fairly representative of the photometric population. 

3. POPULATION SYNTHESIS MODELING 
3.1. Modeling Procedure 

Our goal is to quantify the nebular contribution 
through analysis of the SEDs of a large spectroscopic 
sample at 3.8 < z < 5.0. The advantage of our tech- 
nique is that, for these sources, we can predict the ex- 
act wavelengths of rest-frame optical emission lines and 
thereby remove ambiguities associated with determining 
the photometric redshift simultaneously from contami- 
nated broadband photometry. 

Previous attempts to assess the impact of nebular 
emission on broadband photometry have utilized models 
which include the contributions from both nebular and 
stellar emission. For the reasons outlined in §1, our anal- 
ysis will focus instead on models containing only stellar 
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continuum. The stellar continuum predictions are based 
on the models of Bruzual & Chariot (2003) and the tech- 
nique we will adopt is mostly similar to that described 
in detail in Stark et al (2009). However, we also investi- 
gate how including nebular emission affects the derived 
physical pro perties. To do so , we m ake use of the code 
described in iRobertson et all ()2010D to which the inter- 
ested reader is referred. In this code, line emission is 
calculated from the number of ionizing photons per sec- 
ond, which is provided as output from our population 
synthesis models. 

In the Robertson et al. (2010) code, the intensities 
of hydrogen lines are computed from the values tabu- 
lated in Osterbrock & Ferland (2006), assuming case B 
recombination. We compute the intensities of the lines 
of common metallic species from the empirical results of 
Anders et al. (2003), assuming a gas phase metallicity 
of Z = 0.2 Zq, similar to that measured for galaxies at 
these redshifts (Maiolino et al. 2008, Jones et al. 2012). 
The recombination line luminosities are calculated as- 
suming that the ionising photon escape fraction, f csc , is 
0.2 (see Shapley et al. 2011 for a discussion of expected 
escape fractions). Since we do not rely on the nebular 
models for our primary conclusions, this assumption does 
not affect our results. The continuum contribution from 
bound-free, free-free, and two photon continuum emis- 
sion is also calculated following Osterbrock & Ferland 
(2006). The full nebular template is then added to each 
stellar continuum model, which is then used to calculate 
the synthetic fluxes used in our SED fitting code. 

Since our sample has the virtue of precise spectroscopic 
redshifts, we do not fit the photometric bands spanning 
the Lya forest and Lya emission lines, both of which vary 
significantly from source to source at any given redshift. 
For the redshift range which we are primarily interested 
in (3.8 < z < 5.0), this leaves 7-8 (largely) indepen- 
dent photometric constraints on the SED in GOODS-S 
and 6 constraints on GOODS-N. For consistency with 
the earlier literature, we consider a Salpeter (1955) ini- 
tial mass function with 0.1-100 Mq. Given the rela- 
tively small number of constraints on the SED, we uti- 
lize a moderately restricted grid, varying only the age, 
dust reddening, and normalization factor. We fix the 
star formation history as either constant or rising with 
time following the t power law inferred in Papovich 
et al. (2011). This restricted grid of SFH is supported 
by the results of Reddy et al. (2012) that demonstrate 
that at z ~ 2, the star formation rates inferred from 
exponentially-declining star formation history models do 
not agree with those measured from the observed IR and 
UV fluxes. Nevertheless we have verified that our re- 
sults would not be affected if we had adopted exponen- 
tial decay models. Finally we utilize sub-solar metallic- 
ity (Z=0.2 Zq) motivated by the observations discussed 
above. 

We allow the differential extinction, E(B-V) sta rs, to 
range between 0.00 and 0.50 in steps of 0.02, and we 
limit the model ages to lie between 5 Myr and the age of 
the universe at the redshift of interest. The precise form 
of the dust attenuation curve is, of course, not known 
at z > 4, but we consider the reddening law appropri- 
ate for local starbursts (e.g.,Meurer et al. 1999, Calzetti 
et al. 2000) and a steeper attenuation curve that is ap- 
propriate for the SMC (e.g., Gordon & Clayton 1998). 



The latter appears to be appropriate for young galaxies 
(<100 Myr) at high redshift (Siana et al. 2008, Reddy 
et al. 2010, 2012), a population which might become 
increasingly dominant at z > 4. 

The relative extinction provided to stars and nebular 
emission is not definitively understood at high redshift. 
Expectations from nearby galaxies suggest that the neb- 
ular gas is preferentially more extincted than the stellar 
continuum, as expected if the HII regions lie in dustier re- 
gions than the stars contributing to the integrated stellar 
continuum. Based on observations of local star-forming 
galaxies and starbursts, Calzetti et al. (2000) suggest 
that Ay, neb = Av,sed/0.44. Whether or not this rela- 
tionship holds at high redshift is unclear. Some of the 
first studies of Ha emission in star forming galaxies z ~ 2 
indicated that the nebular gas and stars might be equally 
attenuated (e.g., Erb et al. 2006, Reddy et al. 2010). 
But more recently, new studies have emerged which sup- 
port a factor ~ 2 higher extinction toward HII regions 
(e.g., Forster-Schreiber et al. 2009, Onodera et al. 2010, 
Mancini et al. 2011, Wuyts et al. 2011), similar to that 
observed locally. Clearly an improved understanding of 
how the relative distribution of stars and HII regions de- 
pends on age and mass would greatly benefit attempts to 
simultaneously fit stellar and nebular emission via pop- 
ulation synthesis models. In the nebular+stellar mod- 
els presented in this paper, we will simply assume that 

Ay, SED = Av : nob- 

In the following, we will fit the data with both the 
nebular+stellar and stellar continuum models. For the 
latter, the cleanest method is obtained by fitting the data 
excluding the [3.6] flux measurement, given that this 
band could be contaminated by Ha. However, exclud- 
ing this band means that only one filter is available to 
constrain the SED beyond the Balmer Break. We there- 
fore also fit the data using all photometric information, 
i.e. including the [3.6] measurement. This is discussed 
further in §3.2. For each galaxy, we compute the model 
age, normalization, and Ay which provide acceptable fits 
to the data. The normalization is then mapped to the 
star formation rate and stellar mass appropriate for the 
given template. We compute uncertainties on these pa- 
rameters by bootstrap resampling the data within the 
allowed photometric uncertainties. 

3.2. Nebular Line Strengths 

We infer Ha emission line strengths in our sample of 
3.8 < z < 5.0 galaxies by comparing the observed flux 
density in the [3.6] bandpass to the flux density in that 
filter expected from stellar continuum alone. We explore 
the method which produces the most accurate equivalent 
widths below. 

To verify the reliability of using the broadband flux ex- 
cesses to derive emission line strengths, we examine the 
SEDs of a moderate redshift sample of nebular line emit- 
ters with spectroscopically-measured [OIII] line fluxes 
from WFC3 grism observations of the Hubble Ultra Deep 
Field (Trump et al. 2011). We choose this sample 
rather than larger ground-based samples because use of 
2D grism spectroscopy avoids uncertainties owing to slit 
losses. To characterize the contribution of the emission 
lines to the broadband SEDs, we measure optical through 
mid-IR photometry using the UDF dataset (see McLure 
et al. 2011 for details) and perform population synthesis 
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Fig. 2. — Verification of the flux excess method of measuring emission lines strengths. The four panels show SEDs of spectroscopically- 
confirmed ga l axies with 1.8 < z < 2.3 for which measurements of emission line fluxes are available from the WFC3 grism study of 
ITrump et al.l 1120111) . The flux in the broadband filter contaminated by [OIII] is greater than that expected from the best-fitting stellar 
continuum models in each of the four SEDs. The two displayed stellar continuum models correspond to fits excluding the contaminated 
filter (blue bottom curve) and fits including all available optical and near-IR filters (red top curve). The flux density obtained by subtracting 
the directly-measured emission line contribution to the contaminated broadband flux is shown with a triangle. In each panel, we provide 
the [OIII] emission line flux measured with the WFC3 grism (f spC c) and the the emission line flux inferred from the photometric excess of 
the contaminated filter with respect to the stellar continuum models (fsED) m units of erg/cm 2 /s. 



modelling following exactly the same procedure as de- 
scribed above. [OIII] will either fall in the Ji25-band (at 
1.3 < z < 1.8) or Hieo-band (at 1.8 < z < 2.3). We 
characterize the likely strength of the emission lines by 
comparing the observed broadband flux (in the contami- 
nated filter) to the stellar continuum flux expected from 
the best-fitting population synthesis models. We con- 
sider the stellar models with and without the contami- 
nated bandpass included. In order to ensure a reliable 
measure of the stellar continuum in the vicinity of [OIII] , 
we do not consider galaxies with strong emission lines in 
adjacent infrared filters (e.g., [OIII] in Hi6o and Ha in 
K s ) or those undetected in K s -band. 

We focus our analysis on the four systems in this re- 
maining subset for which the measured [OIII] flux is pre- 
dicted to make a significant (e.g., > 4%) contribution 
to the broadband photometry. In each of these systems, 
the contaminated filter reveals an excess with respect 
to the best- fitting stellar continuum model (Figure 2). 
The fluxes required to produce the broadband excesses 
in the contaminated filter agree well with those mea- 
sured spectroscopically (Figure 2), with the results from 
the two fitting methods typically bracketing the observed 
line flux. While both methods produce remarkably good 
agreement, the line fluxes are slightly more accurate (av- 



erage flux uncertainty of ~ 20%) when the stellar con- 
tinuum is estimated from the fit to the SED including 
the contaminated filter. When the contaminated filter is 
excluded from the fitting process, the inferred line flux is 
typically 1.5-2. Ox greater than measured with the WFC3 
grism. Given the SEDs are fairly poorly sampled in the 
wavelength range where the continuum flux is required, it 
is conceivable that by excluding the contaminated filter, 
the fitting process will prefer redder models which fit the 
flux in the adjacent filters but underpredict the contin- 
uum in vicinity of the emission line of interest. Regard- 
less of the precise reason, it is clear from Fig. 2 that by 
considering the results from both methods, fairly accu- 
rate line strengths can be extracted from the photometry. 

The results of this test therefore motivate use of 
the broadband flux excesses to infer line strengths in 
carefully-selected spectroscopic samples at higher red- 
shifts where direct spectroscopic measurements of neb- 
ular line fluxes are not yet available. As a result of the 
higher redshifts (and the corresponding 1 + z boost in 
observed equivalent width), we expect the nebular line 
contribution to broadband fluxes to be greater than typ- 
ically observed at intermediate redshifts (e.g., Trump 
et al. 2011), allowing the flux excesses to more con- 
sistently stand out with respect to photometric uncer- 



6 



Stark et al. 



35 n 




-1.0 -0.5 0.0 0.5 1.0 
[3.6H4.5] 

Fig. 3. — The effect of emission lines on broadband colors. The 
distribution of [3.6]- [4.5] colours for galaxies at 3.8 < z < 5.0 (blue 
filled histogram) and at 3.1 < z < 3.6 (red shaded histogram). 
The colours in the lower redshift sample reflect the reddened stellar 
continuum, as both filters are free from strong emission lines. In 
contrast, the colors of the 3.8 < z < 5.0 galaxies are shifted toward 
bluer values by Ha emission (which lies in the 3.6^m bandpass). 

tainties. Based on the results in this section, we expect 
the the true line fluxes to be bracketed by our two fitting 
methods, with the most accurate measurements obtained 
by fits to the entire SED. 

4. RESULTS 

We now discuss the results derived from applying our 
technique to the SEDs of spectroscopically-confirmed 
galaxies in the redshift range 3.8 < z < 5.0 over which 
Ha may contaminate the IRAC 3.6/im filter. We com- 
pare the observed [3.6] flux densities to the stellar con- 
tinuum expected from population synthesis models and 
use the results to infer an empirically-based Ha equiv- 
alent width distribution (§4.1). We discuss the possible 
redshift evolution of the nebular line strengths in §4.2. 
Using the empirically-derived equivalent width distribu- 
tions, we examine how nebular emission affects the de- 
rived physical properties of z > 4 galaxies (§4.3). 

4.1. Strength of Nebular Emission Lines 

We begin by comparing the [3.6]- [4.5] colour distribu- 
tion for galaxies at 3.8 < z < 5.0 with that for galaxies 
at 3.1 < z < 3.6 (a redshift range over which the IRAC 
colours are uncontaminated by strong nebular emission). 
This should reveal the impact of Ha emission on broad- 
band fluxes at z > 3. We apply this test for 45 galaxies 
at 3.8 < z < 5.0 with robust flux measurements in the 
4.5/iin filter which constrains the rest-optical stellar con- 
tinuum. The results (Figure 3) point to a significant 
contribution from Ha. The median [3. 6] -[4. 5] colour at 
3.8 < z < 5.0 is 0.33 mag bluer than the median value 
at 3.1 < z < 3.6, consistent with expectations if Ha pol- 
lutes the 3.6/im filter in the higher redshift bin. Note the 



slightly red [3.6]- [4.5] colors of the 3.1 < z < 3.6 sam- 
ple are exactly what is expected for moderately reddened 
(E[B — V] ~ 0.1) galaxies with a constant star forma- 
tion history and luminosity-weighted ages of 100 Myr. A 
Kolmogorov- Smirnoff test demonstrates with confidence 
that these two color distributions are distinct (KS statis- 
tic of D=0.54 with an associated probability by chance of 
8xl0 -8 ). We also consider whether the change in color 
might be due to photometric scatter from increased pho- 
tometric error in the higher redshift bin. We test this 
by randomly perturbing the 3.1 < z < 3.6 [3.6]-[4.5] 
colour distribution according to the IRAC flux errors of 
the 3.8 < z < 5.0 sample. While this can slightly broaden 
the width of the colour distribution, it does not shift the 
median colour to bluer values as observed. 

While the most natural interpretation of the system- 
atic offset is the presence of Ha in the [3.6] filter, it is con- 
ceivable that other effects could contribute. For example, 
one might expect that a systematic offset in [3.6]-[4.5] 
colours might arise from the slightly different rest-frame 
wavelengths sampled and the (potentially) younger ages 
in the higher redshift bin. Examination of population 
synthesis models indicates that intrinsic galaxy evolution 
is not likely to dominate the shift in [3.6]-[4.5] colors. 
Given the median reddening and ages inferred for the 
3.1 < z < 3.6 and 3.8 < z < 5.0 spectroscopic samples, 
we would expect to see [3.6]-[4.5] ~ 0.1. This is simi- 
lar to that observed at 3.1 < z < 3.6, but significantly 
redder than that observed in the 3.8 < z < 5.0 redshift 
range with Ha contamination. We therefore conclude 
that nebular contamination is likely the dominant cause 
of the differences in the [3.6]- [4.5] colours. 

A particularly convincing verification of the above sta- 
tistical test is the fact that we can directly see evidence 
of strong nebular emission in individual SEDs (Figure 
4). Clearly in these examples the flux in the 3.6/im fil- 
ter is not only in excess of that at 4.5/im but is also 
significantly in excess of the stellar continuum of the 
best-fitting population synthesis models. The SEDs of 
galaxies in this redshift range are (not surprisingly) typ- 
ically better fit by models including nebular emission 
(blue lines in bottom panel cf. red lines in top panel 
for stellar continuum models in Figure 4). 

To estimate the strength of Ha, we compute the 
amount by which the observed 3.6/im flux exceeds the 
predicted stellar continuum flux. We define the 3.6/tm 
excess, A[3.6] as the difference between the [3.6] magni- 
tude expected from stellar continuum models which fit 
the SED and the observed [3.6] magnitude. Positive val- 
ues indicate that the observed flux is greater than can be 
accommodated by stellar continuum. This test requires 
an accurate measure of the stellar continuum in the rest- 
optical. Again we limit our sample to those galaxies with 
confident [4.5] detections, as this filter (devoid of strong 
emission lines) is necessary to anchor the population syn- 
thesis models beyond the Balmer Break. 

The distribution of 3.6/im magnitude excesses in our 
spectroscopic sample (Figure 5) reveals that 96% of 
galaxies are observed to be brighter at [3.6] then pre- 
dicted from the best-fitting stellar continuum models. 
The median excess, 0.27 mag, suggests that the typical 
rest-frame emission line equivalent width in the 3.6/^m 
filter at 3.8 < z < 5.0 is 360-450 A. Emission lines there- 
fore contribute nearly 30% of the observed [3.6] broad- 
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Fig. 4. — SEDs of spectroscopically confirmed galaxies at 3.8 < z < 5.0 fit using population synthesis models containing both stellar 
continuum (top row) and stellar+nebular emission (bottom row). Many galaxies in this redshift range show blue [3.6]-[4.5] colours, with 
the 3.6/^m flux significantly in excess of the best-fitting stellar continuum. This flux excess is strongly suggestive of Ha nebular line 
contamination. Not surprisingly, models containing nebular emission provide significantly better fits to the observed photometry, as clearly 
indicated by the agreement between the synthetic (open blue diamonds) and observed (solid black circles) in the bottom row. In the 
following, we will only consider objects for which the Balmer Break can be anchored by a significant (S/N > 5) [4.5] detection, removing 
fainter objects (like that in the right panel) for which an accurate flux excess is difficult to extract. 



band photometry. These values are derived from stel- 
lar continuum model fits that include the contaminated 
3.6/im filter in the modelling. As we demonstrated in §3, 
this method produces the most accurate flux estimates. 
If the contaminated hltcr is excluded from the modelling 
procedure, the inferred continuum level is typically re- 
duced, resulting in a slightly larger median [3.6] excess 
(0.37 mag) and total equivalent widths (520-650 A). Note 
these represent the equivalent width of all emission lines 
in the [3.6] filter. We derive the fractional contribution 
of Ha below. Based on the discussion in §3.2, it is likely 
that these measurements bracket the range of mean neb- 
ular emission line strengths in 3.8 < z < 5.0 galaxies. 
Reassuringly, this EW range is consistent with that re- 
quired to explain the 0.33 mag offset in median [3.6]-[4.5] 
colors of 3.1 < z < 3.6 and 3.8 < z < 5.0 galaxies in Fig- 
ure 3. 

To estimate the equivalent width distribution of emis- 
sion lines contaminating the IRAC 3.6/im filter, we need 
to admit a range of equivalent widths to reproduce the 
observed 3.6/xm photometric excess distribution of Figure 
5a. We assume that equivalent widths are distributed in 
a log-normal fashion, similar to that seen from Ha emis- 
sion locally and at moderate redshifts (e.g., Lee et al. 
2007, Ly et al. 2011). We consider a large grid spanning 
a range of a and /i, the width and mean of the equiv- 
alent width distribution. We translate each equivalent 
width into a 3.6/im excess, applying a photometric scat- 
ter of 20% (a conservative estimate for the average 3.6/im 



magnitude error) and compute the flux excess distribu- 
tion expected from the input equivalent width distribu- 
tion. We find that the observed flux excess distribution is 
well-fit by an equivalent width distribution with a = 0.25 
and < logio(Wr 3 .6]/A) >=2.57 (Figure 5b). This equiv- 
alent width will surely be dominated by Ha emission, 
but other emission lines ( [SII] , [Nil] ) may of course con- 
tribute. The contribution of other lines will depend on 
the physical properties (e.g., metallicity) of the galaxies. 
Our sub-solar (0.2 Z Q ) metallicity models indicate that 
Ha should contribute ~ 76% of the observed equivalent 
width. In this case, the typical Ha EW at 3.8 < z < 5.0 
is < logio(WHa/A) >=2.45. With these assumptions, 
if the [3.6] filter is excluded from the fitting, we find 
< logio(WH Q /A) >=2.61. As above, we adopt this as an 
upper bound to the average Ha EW. 

The level of Ha emission quoted above is actually very 
reasonable given the typical properties of z ~ 4 — 5 LBGs 
(e.g., Stark et al. 2009, Gonzalez et al. 2011a). For 
constant star formation, ionizing photon escape fraction 
of 0.2, and ages of ~ 100 — 250 Myr, we would expect Ha 
EWs to be ~ 200 - 300 A, similar to the range we infer. 
So the observation of [3.6] excesses of 0.2-0.3 mag relative 
to stellar continuum (Figure 3) is exactly what we would 
expect given the shape of the overall SEDs. Indeed, the 
absence of any nebular contamination at 3.8 < z < 5.0 
would have been a far more surprising finding. 

4.2. Evolution of Nebular EW Distribution 
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Fig. 5. — Left: The distribution of 3.6/^m magnitude excesses (A[3.6]) in our 3.8 < z < 5.0 galaxy sample. The magnitude excess is 
defined as the difference between the [3.6] magnitude inferred from the stellar continuum of the best-fitting population synthesis model and 
the [3.6] magnitude observed with Spitzer /IRAC. The positive A [3. 6] values exhibited by our sample indicate that the stellar continuum is 
unable to account for the observed flux in the 3.6 (im filter. The magnitude excess distribution is well fit by a log- normal equivalent width 
distribution with < logioW >=2.57 and a = 0.25 (red curve). Right: Distribution of equivalent widths required to reproduce the observed 
flux excesses (blue histogram) compared to the functional form we adopt for the equivalent width distribution (red curve). 



Before evaluating how nebular emission affects the de- 
rived stellar masses at 4 < z < 7, it is interesting to con- 
sider whether the nebular equivalent width distribution 
is likely to evolve with redshift. Fumagalli et al. (2012) 
recently examined the evolution of Ha equivalent widths 
at lower redshifts, finding that the evolution could be fit 
by a power law oc (1 + z) 18 for galaxies with 10 10 Mq 
in stellar mass (Figure 6). While our sample size is too 
modest to permit a detailed comparison with this trend 
at z > 4, it is of interest to consider how the EWs we 
derived in §4.1 compare to those at lower redshift. We 
note that the Ha equivalent widths presented in Fuma- 
galli et al. (2012) include the contribution of [Nil], while 
the nebular lines strengths we infer from photometric ex- 
cesses include the contribution of all emission lines con- 
taminating the [3.6] filter. Using our 0.2 Z Q population 
synthesis models, we estimate that 82% of the equiva- 
lent width inferred from [3.6] photometric excesses arises 
from Ha and [Nil]. Note that this is slightly larger than 
the percentage estimated in the previous section owing 
to the addition of [Nil] to the calculation. 

Applying this factor to the mean equivalent widths 
presented in §4.1, we compare the Ha+[NII] equivalent 
widths at 3.8 < z < 5.0 to those at z < 2 (Figure 6). It 
is clear that the line strengths derived at 3.8 < z < 5.0 
are consistent with a continued increase in the Ha EW 
in the range 2 < z < 5. While determination of the ex- 
act rate of increase is beyond the scope of this work, we 
note that the Ha EWs at 3.8 < z < 5.0 lie in the range 
expected by simple extrapolations of the power laws de- 
rived in Fumagalli et al. (2012). This increase in Ha EW 
over 2 < z < 5, albeit tentative in nature, is supportive 
of an increase in the sSFR over z > 2, consistent with 
more recent derivations (e.g., Bouwens et al. 2012). 

Given these results, it is reasonable to expect nebular 
lines to be even stronger at z ~ 6 — 7. Unfortunately, 
with both IRAC filters contaminated by strong emission 



lines at these redshifts (Figure 1), we do not have a direct 
method of estimating the nebular line contamination in 
this regime. We thus will consider two cases in the fol- 
lowing sections. First, we assume conservatively that the 
nebular line strengths remain fixed at the values derived 
at 3.8 < z < 5.0. While a fixed EW might seem unlikely 
in light of the power law evolution at lower redshifts, we 
note that the rate of increase in the EW might slow if 
star formation histories transition into a phase of rapidly 
rising star formation rates (e.g., Finlator et al. 2011) at 
5 < z < 7. As a modest upper limit, we consider the 
case whereby the nebular equivalent widths continue to 
increase following a (1+z) 18 power law. We conserva- 
tively adopt the mean EW of the fitting method includ- 
ing the contaminated [3.6] filter as the 3.8 < z < 5.0 
reference value for this upper bound. 

4.3. Effect on z > 4 Stellar Masses 

In §4.1, we demonstrated that strong nebular line emis- 
sion lines make a significant contribution to the broad- 
band flux measurements at z > 3. If these lines are 
not accounted for in population synthesis modelling, the 
rest-optical stellar continuum (and thus the inferred stel- 
lar mass and age) will clearly be overestimated (e.g., 
Schaerer & de Barros 2010). In principle, these issues can 
be addressed through nebular +stellar population synthe- 
sis models described in §3. The drawback of this ap- 
proach is that the 'appropriate' flux from nebular emis- 
sion for any given model is very uncertain, depending not 
only on the escape fraction of ionizing radiation but also 
on the reddening law for the nebular gas, both of which 
are not known at z > 3. 

Here we have attempted to account for these shortcom- 
ings using a method that relies on our empirically-derived 
nebular line equivalent width distribution. For each SED 
we wish to fit, we draw a large number (TV ~ 10 4 ) of Ha 
emission line equivalent widths from the distribution we 
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Fig. 6. — Evolution of the mean Ha+[NII] EW with redshift. 
The values at z < 4 are as compiled in Fumagalli et al. (2012) for 
star forming galaxies with stellar mass in the range logioM*=10.0- 
10.5. The dashed lines show the power law that Fumagalli et al 
(2012) fit to the EW evolution at z < 2 while the dotted lines 
show an extrapolation of this power law to z > 2. The open sym- 
bols denote the EWs appropriate for entire galaxy population in 
Fumagalli et al. (2012), while the solid symbols denote estimates 
for the star-forming subset. The red circles show EWs inferred 
from the photometric excess method in this paper, with the lower 
value arising from SED fits including all filters, and the upper value 
derived from fits excluding the contaminated [3.6] filter (see §3.2 
for details). The range of EWs illustrated at z ~ 6 — 7 repre- 
sents nebular line strengths that we apply to SEDs in §4.3, with 
the lower limit assuming that nebular line strengths remain fixed 
with increasing redshift and the upper limit assuming they follow 
a (1+z) 1,8 power law. 

derived in §4.1 (e.g., Figure 5b). The contribution from 
[OIII] and H/3 (which contaminate the Spitzer/IRAC fil- 
ters at z > 5) is obtained by scaling the Ha EW by 1.7- 
2.0 x (the exact value chosen at random from a uniform 
distribution), consistent with the flux ratios observed in 
sub-solar galaxies at z ~ 2 — 3 (e.g., Hainline et al. 2009, 
Bian et al. 2010, Erb et al. 2010, Richard et al. 2011) 
and those predicted by our nebular+stellar population 
synthesis models. If the galaxy's redshift places any of 
these strong lines in the broadband filters we are fitting, 
we subtract the predicted nebular flux from the photom- 
etry. 

In order to evaluate how the average properties of the 
various z ~ 4 — 7 dropout populations arc affected by 
nebular emission, we consider composite SEDs for the 
B, V, i', and z-band dropout populations (e.g., Stark et 
al. 2009, Labbe et al. 2010b, Gonzalez et al. 2011b) 
binned by rest-UV magnitude. For the purposes of this 
section, we limit our analysis to the most recent deter- 
minations of the composite SEDs (Labbe et al. 2010b, 
Gonzalez et al. 2011b), both of which take advantage 
of WFC3 photometry in the UDF and GOODS fields. 
Since the effect of nebular emission on derived physical 
properties is clearly redshift-dependent (c.f. Figure 1), 
we must account for the distribution of redshifts within 
a particular dropout sample. Thus for each realization 
of the nebular line EW distribution, we also select a red- 
shift from the expected photometric redshift distribution 
of the dropout population under consideration (see e.g., 



Bouwens et al. 2012). We fit these realizations of the 
composite SEDs using the stellar continuum models de- 
scribed in §3. Physical properties are determined in a 
similar manner to that described in §3, with uncertain- 
ties derived from the la intervals of the large number of 
realizations of the equivalent width distribution. 

The impact of nebular emission is clearly seen in the log 
M* - Muv scaling relations shown in Figure 7. These are 
determined for each dropout population with and with- 
out the nebular correction from the composite SEDs dis- 
cussed above. The absolute magnitudes are unchanged 
in this analysis, so the changes shown are due only to 
the contamination of broadband light by nebular emis- 
sion lines. As expected, the impact of nebular emission 
is strongest in the range 5 < z < 7 where nebular lines 
contaminate both IRAC filters. In contrast, the nebular 
correction is less severe for the B-drop (z ~ 4) popu- 
lation, since the Spitzer/IRAC 4.5/im filter is devoid of 
strong emission lines throughout the redshift range cov- 
ered by B-drops (Figure 1). 

Assuming the nebular line EW distribution at z > 5 
remains identical to that determined at 3.8 < z < 5.0, we 
find that the average stellar masses are reduced by factors 
of xl.l, 1.3, 1.6, and 2.4 for the dropout populations 
centered at z ~ 4, 5, 6, and 7 respectively. If the nebular 
line strengths increase with redshift at z > 5 following a 
(1+z) 18 power law (Figure 6; Fumagalli et al. 2012), the 
typical stellar masses are reduced by xl.9 and 4.4 at z ~ 
6 and 7 respectively. We emphasize that these represent 
average corrections applicable to the B, V, i , and z-band 
dropout populations. Individual galaxies throughout this 
redshift range will of course be affected differently. 

To summarize, we have used the equivalent width dis- 
tributions derived in §4.1 to compute the likely contribu- 
tion of nebular emission to broadband photometry. We 
find that stellar masses at z > 6 need to be revised down- 
ward by 2-4 x, with the precise correction depending on 
whether the EW of Ha and [OIII] emission continue to 
increase with redshift beyond z ~ 5. This result has an 
important effect on the log M* - Muv scaling relation, 
which previously was thought to be largely constant with 
redshift at z > 4 (e.g., Stark et al. 2009, Gonzalez et al. 
2011, McLure et al. 2011). It is actually clearly evident 
in the uncorrected log M* - Muv relations presented in 
Figure 7 that without nebular corrections, the M*/Luv 
ratios implied by the composite SEDs increase with red- 
shift. Our analysis indicates that this finding is likely to 
be an artifact of nebular contamination. After correcting 
for line emission, we demonstrate that the M^/Luv ra- 
tios are likely to decrease by x 1.4-2.5 with redshift over 
4 < z < 7. This result has important implications for 
derivation of the stellar mass density (§5.1) and specific 
star formation rate evolution (§5.2) discussed below. 

4.4. Effect of Nebular Continuum Emission 

Nebular continuum emission can also contribute to the 
observed broadband flux density. Most importantly, the 
addition of nebular continuum reddens the intrinsic spec- 
trum at young ages, thereby requiring less dust extinc- 
tion to reproduce the observed colors. This can lead to 
a reduction in the derived star formation rates, in ad- 
dition to the reduction in the stellar masses discussed 
earlier. The effects of nebular continuum can be seen in 
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Fig. 7. — The impact of nebular emission on stellar mass derived using empirical determination of Ha EW distribution at 3.8 < z < 5.0. 
Solid circles show the log M* - Mtjv,1500 relationship corrected for nebular contamination of broadband fluxes following the procedure 
discussed in §4.2. The dashed line gives the best linear fit to these data points. Solid triangles show the relationship with no correction for 
nebular emission assuming that the stellar continuum dominates the broadband flux. Without accounting for nebular emission, the log M 4 
- Mtjv,i500 relationship does not evolve much with redshift. However, incorporating nebular corrections, the normalization decreases by 
1.4-2.5X over 4 < z < 7. The open squares correspond to nebular corrections derived assuming that the EW of nebular emission increases 
with redshift following the power law shown in Figure 6. 



the nebular+stellar models, which typically display spec- 
tral edges in the vicinity of the Balmer/4000 A break 
(bottom panel of Figure 4). Whether significant neb- 
ular continuum is actually contributing to the spec- 
tra of high-z galaxies is unclear. Without better sam- 
pled SEDs (from e.g., medium-band near-IR photome- 
try) that might probe such spectral edges, it is difficult 
to verify the presence of nebular continuum, as we are 
able to with nebular emission lines. Nevertheless if emis- 
sion lines are very strong, it is likely that there is also 
a significant contribution from nebular continuum emis- 
sion. 



To quantify the likely impact of nebular continuum on 
the derived physical properties, we examine how the in- 
ferred dust reddening (and SFR) are affected when the 
nebular continuum is added. To do so, we compare the 
dust attenuation necessary to reproduce a fixed UV con- 
tinuum slope for models with and without nebular con- 
tinuum emission included. Assuming a Calzetti redden- 
ing law, we calculate the the dust attenuation that repro- 
duces a UV slope with (3 = —1.5 as a function of model 
age. Owing to the redder intrinsic slopes, the inferred 
dust attenuation is reduced for the nebular+stellar mod- 
els. The effect is most pronounced at the youngest ages, 
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Fig. 8. — Evolution in the Stellar Mass Density: The mass density is computed by integrating the stellar mass function to a fixed UV 
luminosity limit (left panel) and fixed stellar mass limit (right panel). The main advance with respect to earlier work is the inclusion of 
corrections for nebular emission contamination of Spitzer/IRAC filters. These are computed using the nebular EW distribution derived 
from our spectroscopic sample in §4.1 and assume the Ha EW distribution continues to evolve as a power law at z > 5. The stellar mass 
functions are derived by combining the nebular-corrected log M*-Ltjv relation with the UV luminosity functions presented in Bouwens 
et al. 2012. We assume that the scatter in the log M*-Ltjv relation is a ~ 0.5, consistent with previous studies (e.g., Gonzalez et al. 
2011a). The blue swath in the right panel shows the stellar mass density implied by the evolving star formation rate density, computed by 
integrating the UV luminosity functions (see Robertson et al. 2010 for details). 



with the inferred attenuation up to 20% lower for nebu- 
lar+stellar models with ages <30 Myr. As a result of the 
reduced attenuation, a lower normalization is required to 
match the observed flux density, bringing down the in- 
ferred star formation rates by up to the ~ 20% level for 
the youngest systems. In practice, the impact of neb- 
ular continuum is more complicated and non-trivial to 
predict, depending strongly on the shape (i.e. age) of 
the observed SED. Given these uncertainties, our anal- 
ysis in the following sections will focus on how physical 
properties are affected by nebular emission lines. 

5. DISCUSSION 

In the previous section, we used the broadband SEDs 
of a large sample of spectroscopically-confirmed galaxies 
to infer the distribution of nebular line strengths in UV- 
selected galaxies at 3.8 < z < 5.0. We showed that the 
stellar masses inferred from population synthesis mod- 
elling are reduced at z > 5 when the contribution of 
these lines to broadband flux densities is removed. In 
this section, we consider the implications of these results 
for our current picture of early mass assembly. 

5.1. Stellar Mass Density at z > 3 

In §4.2, we quantified the extent to which the stellar 
masses of z > 3 galaxies are affected by nebular emission. 
Here, we seek to utilize these results to estimate the stel- 
lar mass density (SMD) evolution at z > 3. To derive 
the mass densities, we combine the log M* - Mtjv rela- 
tionship with UV luminosity functions (LFs) in a manner 
mostly similar to that outlined in Gonzalez et al. (2011). 
Briefly, we extract a large number (N ~ 10 5 ) of lumi- 
nosities from the measured UV LFs (e.g., Bouwens et al. 
2011). We convert these luminosities to stellar masses 
using the log - Mtjv relationship, and an estimate of 



the scatter about the median. Whereas earlier studies 
held the log M* - Mtjv relationship fixed with redshift at 
z > 4, the strong redshift dependence of nebular contam- 
ination (Figure 1) forces us to reconsider the evolution 
of M^/Ltjv ratios with redshift. 

We compute the slope and normalization of the z ~ 4 
log M* - Mtjv relationship using the large sample of 
LBGs discussed in Stark et al. (2009). For simplicity, we 
assume the slope remains constant at z > 4 and consider 
only evolution in the normalization of the relationship. 
To compute the zero-points of the log M* - Mtjv relation 
at z ~ 5, 6, and 7, we adjust the measured z ~ 4 relation 
to account for the relative normalisation of the nebular- 
corrected log M* - Mtjv relationships shown in Figure 7. 
To obtain a tentative estimate of the z ~ 8 stellar mass 
density, we apply the z ~ 7 log M* - Mtjv relationship 
to the z ~ 8 UV LF. In all cases, we use the nebular 
corrections derived assuming an evolving Ha EW distri- 
bution (see Figure 7), but we also discuss how these re- 
sults would change if the EW distribution remains fixed 
at z > 5. 

In addition to measurement of the log M* - Mtjv rela- 
tionship, accurate determinations of the dispersion are 
necessary to account for low luminosity galaxies with 
large M* /Ltjv ratios. If scatter is not accounted for, 
the mass functions will be incomplete and mass densities 
(above a fixed mass limit) will be underestimated. While 
a measurement of the observed scatter at z ~ 4 (0.5 dex) 
was made in Gonzalez et al. (2011a), the intrinsic scat- 
ter is likely lower due to systematic uncertainties in the 
modelling as well as the effects of nebular contamination. 
In the following, we assume the intrinsic scatter is in the 
range 0.2-0.5 dex. 

We focus first on the UV luminosity limited measure of 
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the SMD, considering only those galaxies with luminosi- 
ties greater than Muv = — 18. The results are presented 
in Figure 8. As discussed in earlier sections, if nebular 
emission is not accounted for in the modelling, the in- 
ferred M*/Luv ratios actually increase moderately with 
redshift, leading to artificially high stellar mass densities 
at z > 5. The inclusion of nebular emission reduces the 
mass density at z ~ 7 by up to 4 x , while having little 
effect at z ~ 4. As a result of these redshift-dependent 
corrections, the evolution in the SMD is fit by a steeper 
power law ([1+z] -4,7 ) than reported previously. 

We also consider the SMD in galaxies with stellar 
masses in excess of 10 s M©. For consistency with earlier 
measurements of the SMD, the measurements we present 
in Figure 8 assume that the log M* - Muv relationship 
has an intrinsic dispersion of 0.5 dex. If the scatter is 
instead only 0.2 dex, for example, we would find stellar 
mass densities that are ~ 1.6 — 2. Ox lower. 

While previous estimates of the 6 < z < 7 SMD ap- 
peared broadly consistent with the integral of the z > 7 
SFRD (after an appropriate correction for stellar mass 
loss and recycling), the SMD appeared to be on the high 
end of the range implied by the SFRD (e.g., Robertson et 
al. 2010). We compare our revised SMD to those implied 
by the SFRD in the right panel of Figure 8. The mass 
density implied by the SFRD is calculated in a similar 
manner as Robertson et al. (2010), updated to include 
the latest measurements of the UVLF (Bouwens et al. 
2011). With appropriate corrections for nebular emis- 
sion, the mass densities now appear in better agreement 
with the integrated SFRD. 

The stellar mass density provides a useful integral con- 
straint on earlier star formation. Spitzer observations of 
galaxies at z ~ 6 — 7 therefore offer a valuable measure 
of the likely ionizing output of galaxies at 7 < z < 15. 
Therefore as measurements of the stellar mass density 
become more reliable, they will offer insight into the con- 
tribution of galaxies to reionization, complementing in- 
ferences from the UV luminosity function. 

We have demonstrated in this paper that corrections 
for nebular emission are a crucial aspect of obtaining a re- 
liable census of stellar mass in the early universe. In light 
of the reduced mass density required by nebular contam- 
ination, we re-consider the ability of galaxies to achieve 
reionization by z > 6, updating the calculation presented 
in Robertson et al. (2010). Given the consistency with 
the SFRD (Figure 8), these results are not surprisingly 
similar to inferences obtained from the UV LF (e.g ., 
iBouwens et~alll2012at Ikuhlen fc Faucher-Giguerell2012h . 
The UV output implied by the mass density is in prin- 
ciple sufficient to achieve reionization by z ~ 6 — 8 but 
struggles to account for the optical depth to electron scat- 
tering implied by WMAP (e.g., Larson et al. 2011). 

Understanding this photon shortfall will require im- 
proved knowledge of how much star formation occurs 
beyond z ~ 10. While direct detection of z > 10 galax- 
ies will likely have to wait until JWST, Spitzer offers a 
unique means of progress in the coming years. By ob- 
taining stellar mass estimates for the emerging samples 
of z ~ 9 — 10 galaxies (e.g Bouwens et al. 2011b, Zheng 
et al. 2012), it will be possible to obtain some of the first 
constraints on the contribution of galaxies to the cosmic 
ionisation history beyond z ~ 10. 



5.2. sSFR Evolution 

The reduced stellar masses we infer in §4 clearly will 
affect the evolution of the sSFR at z > 4. To estimate the 
impact, we compute the sSFR in fixed stellar mass bins 
using a similar approach as for the stellar mass function. 
We draw a large number (N ~ 10 5 ) of luminosities from 
the latest measures of the z ~ 4 — 7 UV LFs (Bouwens 
et al. 2012a). For each luminosity and redshift bin, we 
compute a stellar mass using the log M*-Muv relation- 
ship derived in §4.1. We consider the case in which the 
strength of nebular emission is constant at z > 4 and 
also that in which the emission line equivalent widths 
increase with redshift (e.g., Figure 6). 

The SFR is computed from Muv through a series of 
steps. We account for dust extinction using the UV con- 
tinuum slopes. For each realisation of the UV LF, we 
draw a UV slope, /3, by adopting the redshift-dependent 
f3 — Muv scaling relationships (Bouwens et al. 2012b). 
The UV slope is then converted to a dust attenuation fac- 
tor at 1600 A via the Meurer et al. (1999) IRX-/3 relation 
(Ai 600 =4. 43+1. 99/3). The UV luminosity is converted to 
SFR following the canonical Madau et al. (1998) and 
Kennicutt et al. (1998) relation L uv = (SFR/M© yr^ 1 ) 
8.0xl0 27 ergs s _1 Hz -1 . This relationship assumes a 0.1- 
125 M Q Salpeter IMF and constant star formation rate 
of > 100 Myr. Finally, by examining the SFR and M* of 
these realizations, we compute the median sSFR of the 
four dropout samples with stellar mass of 5 x 10 9 M Q . Be- 
fore examining the results of this calculation, we discuss 
two important issues that we have hitherto neglected. 

First we consider how the sSFR is affected by scatter 
in the log M*-Muv relationship. Note that the sSFR 
will be overestimated if one merely uses the log M*-Muv 
relation without taking into account the abundant pop- 
ulation of lower SFR objects with large M*/Luv ratios. 
This issue is dealt with in detail in Reddy et al. (2012) 
for galaxies at z ~ 2 — 3. Unfortunately, as we discussed 
in §5.1, the intrinsic scatter is very poorly constrained 
in UV-selected samples at z > 4. As a result, previ- 
ous estimates of the sSFR at z > 4 have not accounted 
for M*/Luv scatter. To estimate how this shortcom- 
ing would affect the sSFR, we add log M^-Muv scatter 
to the LF realization method described above. If the 
0.5 dex observed z ~ 4 scatter reported in Gonzalez et 
al. (2011a) is entirely intrinsic, then the median sSFR 
would be reduced by 2.8 x at z ~ 4. Note that one might 
find slightly different adjustments for the same scatter at 
other rcdshifts owing to evolution in the luminosity func- 
tion. In §5.1, we suggested that the intrinsic scatter is 
likely lower as systematic uncertainties in modelling (in- 
cluding uncertainties in the nebular corrections) surely 
broaden the dispersion in the stellar masses at fixed UV 
luminosity. In this case, fewer low SFR galaxies con- 
tribute to the sSFR distribution at fixed mass, resulting 
in a larger median sSFR. For example, a scatter of 0.2 
dex would translate into a reduction of just 1.2 x with 
respect to the case of no scatter. Physically, one may ex- 
pect that the scatter at fixed luminosity would increase 
somewhat between z ~ 7 and z ~ 4, as galaxies have 
had more time to undergo punctuated episodes of star 
formation which elevate both their star formation rates 
and mass off of the main sequence. We consider these 
possibilities in our discussion below. 
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Fig. 9. — Left: Evolution in the specific star formation rate (sSFR). Our new measurements include stellar masses corrected for 
nebular emission line contamination. The solid red circles at z > 4 show values derived assuming the nebular line EW distribution at 
4 < 2 < 7 remains identical to that derived at 3.8 < z < 5.0 (Figure 5). The solid squares correspond to values obtained when an 
evolving nebular line EW distribution is adopted. The error bars show the assumed scatter about the mean sSFR taken from Reddy et al. 
(2012). Right: Comparison of observed sSFR to contemporary theoretical models. The solid lines show the sSFR evolution predicted from 
cosmological simulations discussed in Dave et al. (2011), with the blue line corresponding to their "slow wind" model and the light green 
line corresponding to the momentum driven wind "vzw" model. These models provide adequate fits at the highest redshifts (z > 5) but 
undershoot the observed values at z ~ 2 — 4. 



We now examine a second issue which pushes the sSFR 
in the opposite direction. In particular, we examine how 
scatter (and perhaps systematic offsets) in the conversion 
between dust-corrected Ltjv and SFR affect our sSFR de- 
termination. The conversion between UV luminosity and 
SFR that we use is valid for galaxies with model ages in 
excess of 100 Myr. Above this age, the conversion fac- 
tor changes little for an assumed constant star formation 
history. But below 100 Myr, a larger SFR is required 
to produce a fixed Ltjv (see Figure 25 of Reddy et al. 
2012). For example, a galaxy with model age of 10 Myr 
requires a 1.8 x larger SFR to reproduce the same Ltjv as 
a galaxy with 100 Myr. With the reduced ages implied by 
the nebular corrections, it seems likely that such young 
systems are present in z > 4 dropout samples. Inclusion 
of dispersion in the model ages will preferentially shift 
the median SFR to larger values, resulting in somewhat 
larger sSFR. Furthermore, it is of course conceivable, if 
not likely, that such young systems will become more 
common both at higher redshift, requiring a systematic 
shift toward higher sSFR at earlier times. 

We now turn to the derived sSFR evolution, which is 
shown in Figure 9a. First, ignoring the effect of neb- 
ular emission and the scatter discussed above, we find 
that the sSFR actually decreases marginally with red- 
shift over 4 < z < 7, similar to the findings of Bouwens 
et al. (2012b). This is driven largely by the redshift- 
dependence of the UV continuum slope (3 versus Mtjv 
relationship. Galaxies at higher redshifts have bluer UV 
continuua (e.g., Bouwens et al. 2012b, Finkelstein et 
al. 2012), reducing the dust-corrected SFR for a given 
Mtjv- Considering the fixed M*/Ltjv ratios assumed in 
previous studies (e.g., Stark et al. 2009, Gonzalez et 
al. 2010), it is straightforward to understand this result. 
As we have discussed above, without nebular corrections, 
the data actually support a mild increase in the M*/Ltjv 



ratios with increasing redshift at z > 4 (Figure 7); if the 
M*/Ltjv ratios weren't held fixed (and nebular emission 
not considered), one would have derived a more rapid 
decrease in sSFR at z > 4. 

Incorporating our corrections for nebular emission re- 
duces the M*/Ltjv ratios in the z ~ 5 — 7 LBG samples, 
increasing the sSFR in this redshift range. If the nebular 
line EW distribution at z > 5 is similar to that seen in 
Figure 5b, we find that the sSFR of galaxies with fixed 
stellar mass begins to show evidence for positive evolu- 
tion with redshift, with the z ~ 7 value (9.6 Gyr -1 ) 4x 
larger than that at z ~ 2. This can be viewed as a con- 
servative estimate of the sSFR evolution. As we have ar- 
gued, however, it more likely that the equivalent width 
of Ha and [OIII] increase in strength with redshift at 
z > 4, consistent with the evolution seen at intermediate 
redshift (Fumagalli et al. 2012). Under these assump- 
tions, the derived sSFR shows greater redshift evolution, 
with the z ~ 7 sSFR (14 Gyr -1 ) roughly 6 x greater 
than that at z ~ 2 (Reddy et al. 2012). While intrinsic 
scatter in the M*/Ltjv ratios might bring these numbers 
down somewhat (perhaps explaining the excess seen at 
z ~ 4), this is likely counteracted somewhat by scatter 
and/or systematic evolution in the SFR/Ltjv ratios and 
possibly a shift toward reduced scatter in the M*/Ltjv 
ratios at higher redshifts. 

To summarize, with the new dust corrections (Bouwens 
et al. 2012b) and adjustments for nebular emission con- 
tamination, we now find evidence for a power law in- 
crease in the sSFR at z > 2 that is much more consistent 
with theoretical expectations than previous observations 
indicated. Both the absolute values and rate of increase 
of the sSFR we derive at z > 5 are very similar to those 
predicted in the simulations of Dave et al. (2011a). In- 
triguingly the sSFR at 2 < z < 4 still remains moder- 
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ately in excess of theoretical expectations. As we have 
discussed above, the z ~ 4 estimate of the sSFR might 
come down somewhat owing to scatter in the M*/Luv- 
But the z ~ 2 — 3 sSFR measurements include the effects 
of scatter, so the discrepancy remains puzzling, especially 
in light of the emerging agreement at z > 5. Previous 
theoretical studies have focused on a variety of explana- 
tions for the tension at z ~ 2, including a time- varying 
initial mass function (e.g., Dave 2008, Narayanan & Dave 
2012). Continued efforts along these lines are required to 
simultaneously explain the high sSFR at z ~ 2 along with 
the revised z > 4 sSFR estimates (Fig. 9). 

Recall that previous indications of a nearly fiat sSFR 
in fixed stellar mass bins at z > 2 required a mecha- 
nism by which star formation is made increasingly ineffi- 
cient at earlier times. Possibilities included the inefficient 
formation of molecular hydrogen in metal-poor galaxies 
(e.g., Robertson&Kravtsov 2008; Gnedin et al. 2009; 
Krumholz & Dekel 2012), or an increase in the mass out- 
flow rate per unit star formation with redshift. The up- 
dated estimates of the z > 3 sSFR no longer obviously re- 
quire a significant suppression of star formation in galax- 
ies with stellar mass in excess of 10 9 M Q . The current 
measurements seem consistent with a picture whereby 
the rapidly increasing baryon accretion rates translate 
into higher sSFR at earlier times. 

As with the stellar mass density, there is room for con- 
siderable improvements to these estimates in the coming 
years. By providing more individual detections of UV- 
faint galaxies, deeper Spitzer data will enable improved 
measurements of the slope and scatter of the log M*-Muv 
relationship in the redshift range considered in Figure 9. 
It is also of interest to extend these measurements to 
z ~ 8. In this redshift regime, [OIII] lies in the 4.5^im 
filter, while the 3.6/im filter is devoid of strong lines. 
Thus, with deep Spitzer data, the SEDs of z ~ 8 systems 
enable a unique method of deciphering how the strength 
of nebular emission evolves over 5 < z < 8, one of the 
key uncertainties in the current analysis. 

6. SUMMARY AND CONCLUSIONS 

Measurement of the evolving stellar mass and sSFR 
distributions have proven critical to our understanding 
of early galaxy assembly and the UV photon budget of 
rcionization-era galaxies. Recently, it has become clear 
that many of these early estimates might be significantly 
in error due to the contamination of the Spitzer/IRAC 
bandpasses by nebular emission lines. Knowledge of the 
strength of these emission lines is necessary for robust de- 
terminations of the stellar mass density and sSFR evo- 
lution. As these emission lines are shifted out of the 
observed atmospheric window, direct spectroscopic mea- 
surements will not be feasible until JWST. 

In this paper, we present a method which enables con- 
straints on nebular emission at z > 4 by combining 
large spectroscopic samples and deep Spitzer photometry. 
Like Shim et al. (2011), we focus on the redshift range 
3.8 < z < 5.0, over which the IRAC [3.6] filter is contami- 
nated by strong emission lines (Ha, [Nil], [SII] ) while the 
[4.5] filter is free of nebular contamination. Examining a 
carefully-selected subset of 45 galaxies, we find that the 
3.6/im flux is systematically in excess of the expected 
stellar continuum flux, revealing the presence of strong 



nebular emission. No excess is seen in a spectroscopic 
sample at 3.1 < z < 3.6, a redshift range over which no 
strong emission lines contaminate the IRAC filters. We 
use the photometric excesses in the contaminated [3.6] 
filter to estimate the equivalent width distribution of Ha 
emission at 3.8 < z < 5.0. Equipped with this measure 
of nebular emission at high-redshift, we re-evaluate the 
evolution in the sSFR and stellar mass density at z > 4. 
Our primary conclusions from this analysis are summa- 
rized below. 

1 . We find that the mean equivalent width of emission 
lines contaminating the [3.6] filter is < logio(W3.6/A) > 
~ 2.57 - 2.73. We estimate that ~ 76% of this signal 
arises from Ha, implying an average Ha equivalent width 
of < logio(W HQ /A) > ~ 2.45 - 2.61 at 3.8 < z < 5.0. 

2. The mean Ha equivalent width inferred at 3.8 < z < 
5.0 appears greater than that for similar star-forming 
samples at lower rcdshifts. While definitive knowledge 
of the evolution in the Ha EW surely awaits direct 
spectroscopic measurement, the evolution we infer over 
2 < z < 5 is certainly consistent with the (1+z) 18 power 
law derived in Fumagalli et al. (2012). This likely re- 
flects an increase in the sSFR at z > 2, and importantly 
suggests that the EW of nebular emission continues to 
increase at z > 5. 

3. Using the Ha EW distribution we derive at 3.8 < 
z < 5.0, we explore how nebular contamination is likely 
to affect the physical properties of galaxies at z > 3. 
We find that the stellar masses are reduced, on aver- 
age, by 1.1, 1.3, 1.6, and 2.4x for dropout samples with 
mean redshifts of z ~ 4, 5, 6, and 7, respectively. If the 
equivalent widths of nebular lines continue to increase 
in amplitude at z > 5, we estimate that the reduction in 
the stellar masses are likely to increase to 1.9 and 4.4 x at 
z ~ 6 and 7. We note that these corrections are represen- 
tative for average measures of the dropout populations 
and not individual galaxies. 

4. As the stellar mass density provides a valuable inte- 
grated measure of early star formation, constraints on the 
level of nebular contamination are critical to our knowl- 
edge of the ionizing photon budget of galaxies through- 
out the reionization era. After correcting for nebular 
emission contamination, we find a factor of ~ 2 x reduc- 
tion from previous estimates. The downward revisions 
to the stellar mass density improve consistency with ex- 
pectations from the integrated star formation rate den- 
sity. Extending such nebular-corrected measurements to 
emerging galaxy samples at z ~ 8— 9 will yield an integral 
constraint on the UV photon budget during z ~ 10 — 15. 

5. Whereas previous derivations showed little evolu- 
tion in the sSFR of fixed mass galaxies over 2 < z < 7, we 
demonstrate that after accounting for nebular emission 
and correcting for dust, the sSFR increases by 4-6 x over 
2 < z < 7. The absolute sSFR values inferred at z > 5 
appear largely similar with predictions from simulations. 
While there certainly remains room for improvement (in 
both the data and the modelling), the increase in the 
sSFR at z > 4 seems consistent with a picture whereby 
increasing baryon accretion rates at larger redshift trans- 
late into larger sSFR in galaxies of a fixed stellar mass. 
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TABLE 1 

GOODS BRIGHT 3.8 < Z < 5.0 SPECTROSCOPIC sample 



ID RA(J2000) DEC(J2000) z spoc z 850 A[3.6] 
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